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ABSTRACT 

It is shown that pure exponential discs in spiral galaxies are capable of supporting 
slowly varying discrete global lopsided modes, which can explain the observed fea- 
tures of lopsidedness in the stellar discs. Using linearized fluid dynamical equations 
with the softened self-gravity and pressure of the perturbation as the collective effect, 
we derive self-consistently a quadratic eigenvalue equation for the lopsided perturba- 
tion in the galactic disc. On solving this, we find that the ground-state mode shows the 
observed characteristics of the lopsidedness in a galactic disc, namely the fractional 
Fourier amplitude Ai increases smoothly with the radius. These lopsided patterns 
precess in the disc with a very slow pattern speed with no preferred sense of preces- 
sion. We show that the lopsided modes in the stellar disc are long-lived because of 
a substantial reduction (~ a factor of 10 compared to the local free precession rate) 
in the differential precession. The numerical solution of the equations shows that the 
ground-state lopsided modes are either very slowly precessing stationary normal mode 
oscillations of the disc or growing modes with a slow growth rate depending on the 
relative importance of the collective effect of the self-gravity. N-body simulations are 
performed to test the spontaneous growth of lopsidedness in a pure stellar disc. Both 
approaches are then compared and interpreted in terms of long-lived global m = 1 
instabilities, with almost zero pattern speed. 

Key words: galaxies: kinematics and dynamics — galaxies: evolution — galaxies: 
spiral — galaxies: structure — galaxies: general. 



1 INTRODUCTION 

There is growing observational evidence that large-scale lop- 
sided asymmetry is a common phenomenon in spiral galax- 
ies, as seen in the stellar distribution (e.g., Block et al. 1994, 
Rix & Zaritsky 1995, Bournaud et al. 2005) as well as in the 
distribution of the atomic hydrogen gas (e.g., Baldwin et 
al. 1980, hence BLS, Richter & Sancisi 1994, Haynes et al. 
1998). These observations show that about 30% of the field 
spirals show significant disc lopsidedness. 

Various theoretical explanations have been put forward 
towards explaining the physical origin of the lopsided mass 
distribution in spiral galaxies. Some of these are: the disc 
response to a distorted halo (Jog 1997, 1999), satellite infall 
onto a galaxy (Zaritsky & Rix 1997), and tidal interaction 
and asymmetric gas accretion onto a galactic disc (Bournaud 
et al. 2005) . Cooperation of orbital streams to the lopsided 
pattern was used to drive the lopsided instabilities in disc 
galaxies (Earn & Lynden-Bell 1996). Lovelace et al. (1999) 
found strongly unstable eccentric motions within one disc 
scalelength, and the model with an off-centred disc w.r.t. 
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the halo by Levine & Sparke (1998) also shows lopsidedness 
in the inner regions- but both these do not give lopsidedness 
in the outer regions where it is seen in galaxies. So far some 
theoretical studies have shown that lopsided instabilities ex- 
ist only in counter-rotating stellar discs in which the fraction 
of retrograde stars is large (Hozumi & Fujiwara 1989; Sell- 
wood & Valluri 1997) . Since counter rotation is a rarely seen 
phenomenon in stellar discs (Kuijken, Fisher, & Merrificld 
1996; Kannappan & Fabricant 2001) it leaves a room for the 
search of lopsided instability in normal differentially rotating 
galactic discs. A kinematic origin for lopsidedness in stellar 
disc indicates a short winding-up time scale which is less 
than a Gyr (BLS). N-body simulations show the disc lop- 
sidedness to be long-lived to ~ 3-4 Gyr, however the physical 
reason for this is not understood. It has been shown that in 
field galaxies, the disc lopsidedness is uncorrelated to the 
strength of tidal interaction (Bournaud et al. 2005), or the 
presence of nearby companions (Wilcots & Prescott 2005). 
The latter point is exemplified by the case of M101 which 
is an isolated galaxy and yet is strongly lopsided. Thus, an 
internal origin for the generation of the disc lopsidedness is 
worth investigating. Because it has the potential to provide 
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answer to the question: why are so many isolated galaxies 
lopsided? This provides the motivation for our present work. 

The two main observed characteristics of the disc lop- 
sidedness are that: first, Ai, the fractional Fourier ampli- 
tude of lopsidedness increases smoothly with radius as seen 
in Rix & Zaritsky (1995), Bournaud et al. (2005), and mea- 
sured between 1.5-2.5 disc scalelengths. Second, the phase of 
lopsidedness is nearly constant with radius (Rix & Zaritsky 
1995, Angiras et al. 2006). The latter indicates that the disc 
lopsidedness is a global feature (e.g., Jog 1997). 

Guided by these observed features, in this paper, we 
study self-consistently the behaviour of the global lopsided 
perturbation in a purely exponential galactic disc using the 
fluid dynamical approach. The lopsidedness emerges as a 
global instability arising due to the collective effect of the 
self-gravity of the perturbation in the disc. The resulting 
pattern speed is almost constant with radius and thus avoids 
the winding up problem due to the differential precession. In 
our picture, the lopsidedness manifests itself as a classic case 
of the negative damping phenomenon in a self-gravitating 
system. Note that the modes with m = 1 symmetry have in 
general no inner Lindblad resonance in the disc. So refrac- 
tion and dissipation of acoustic waves is less of a problem 
(Block et al. 1994). Unlike the usual m = 2 spiral pattern, 
there exists no definite pattern speed measurement for the 
rn — 1 lopsided mode in the galactic disc. 

For simplicity and to separate the various different dy- 
namical mechanisms, we treat the lopsidedness in a purely 
exponential stellar galactic disc and disregard the effect of 
halo. We also postpone the treatment of a dissipative com- 
ponent. In a future paper, we will include the effect of halo 
and treat a similar global calculation that is applicable to 
the atomic hydrogen gas located at larger radii. We tackle 
the existence of global m = 1 mode both from analytical 
development, and through numerical simulations. 

We organise the paper in the following way: In section 
2 we formulate the dynamical equations governing the lop- 
sided mode in the disc and provide a numerical scheme for 
solving the matrix eigenvalue equation. The results from 
this linear approach are described in section 3. Section 4 
describes the N-body simulations of a pure stellar disc, and 
the study of the m = 1 modes. The comparison between the 
two approaches are compared and interpreted in section 5. 
Section 6 summarizes our conclusion. 



cylindrical coordinate system (R, ip) can be written as: 



Dv'r 
Dt 

Dv[, 



■2€lv' 



1 dP' R | K 



PS' 1 d(RT,°v' R ) 
~Dt + R dR 



7777 E° OR ' li^" 
1 d& _ 1 dP' v 
R dtp RE dip 



P' 

^R 



R dip 



= 



(1) 
(2) 
(3) 



In the above equations D/Dt = d/dt + Q.d/dp and k is 
the epicyclic frequency in the disc. We consider two sepa- 
rate cases namely a cold disc and a hot disc throughout this 
paper. In the case of a cold disc, the equlibrium state is 
defined by a balance between the differential rotation and 
the softened gravity. Whereas in the equilibrium state, a 
hot disc is supported by the differential rotation and pres- 
sure against the softened gravity. Since we are primarily in- 
terested in the stability properties of the global lopsided 
modes(m = 1) in the disc, we consider all the perturbed 
variables as X'(R,ip,t) = X' \R)e^ v ~ ut) about the equilib- 
rium state. Stability of the mode depends on the sign of the 
imaginary part of ui (Im(a>)). Lopsided modes in the disc 
become unstable, in other words they are self-excited, when 
Im(w) > . They are stable decaying modes when Im(cj)< 0. 
Modes with Im(u) = are in general stationary van Kam- 
pen modes; outside the continuum, these modes are pure 
normal mode oscillation of the whole disc. 

We study here the slowly varying (u> «Q) global lop- 
sided mode using softened self-gravity and the pressure as 
the dominant collective effect in the disc. In effect we aim 
to study here the lopsidedness in a real galactic disc. The 
effect of softened self-gravity was studied earlier on the slow 
modes in Keplerian discs by Tremaine (2001). We first ne- 
glect the collective effects due to the pressure in the stellar 
disc in order to bring out clearly the effect of self-gravity in 
governing the properties of these large-scale global modes in 
the galactic disc. Then we include the effect of velocity dis- 
persion to study how the behaviour of lopsidedness changes 
in the hot disc. Substituting the above form of the perturbed 
variables into eqs.(l-3) and solving for the velocity field in 
the limit ui « Q we obtain: 



v R = — 



[W' R + 2W' V ] 



(4) 



2 DYNAMICS OF GLOBAL LOPSIDED MODE 

We formulate here the dynamical equations governing the 
global lopsided mode in a thin, initially axisymmetric self- 
gravitating disc having an exponential distribution of surface 
density (E°). By assuming hydrostatic equilibrium along the 
vertical direction, we ignore the vertical structure of the disc 
by integrating all physical quantities over z. The disc is dif- 
ferentially rotating with angular speed Q(R) so that the un- 
perturbed velocity field in the plane is given by v = (0, RQ,). 
We write the velocity field, surface density, potential and the 
pressure of the perturbed disc as : vr = v' R , v v — RQ + v' v , 
E = E° + E', $ = $° + $' and P = P° + P' . 
Then the linearized Euler and continuity equations in the 



^ = w+m=^) ^ Wk+2W ^ (5) 

Where 

, = d& 1 dP' R P'^-P'r 
R dR E° dR i?E° 

w = ^ + 

v R + i?E° 

In the above expressions lu p — — k is the free preces- 
sion frequency of the m = 1 lopsided mode and 8 — k/Q 
in the disc under consideration. Interestingly, note that the 
ratio of the two velocities is no longer a simple scalar as it 
is in the Keplerian limit (5 = 1). 

Substituting the perturbed variables in continuity equation 
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and simplifying the equation by using the limit to << fl we 
get: 
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Now connecting the velocity field (eqs. [4-5]) with the 
continuity eq.(6) and performing straightforward mathemat- 
ical manipulations we arrive at the following compact equa- 
tion: 

oj 2 S' + S(S',$'V + <5(S',$') = (8) 

Where 

£>(£', $') = - [2w p E' + (0\& + OpS')] (9a) 
<S(E', $') = c^E' + uj p (O g & + OlT,') 

-e(R)(0 2 $' + 2 p Z') (96) 

In deriving the above eq.(8) we have used the fact that the 
vertically integrated pressure is a function of the surface den- 
sity E. In particular we have used the perturbed pressures as 
Pr = a R^' an d P!p = o^E'. So we take into account of the 
anisotropy in pressures in the plane of the disc. In eq.(9b) 
e(R) = adWp/dR where a is defined below in eq.(lla). Note 
that e would be zero if there were no differential precession 
in the disc. 

We call 0\ and O 2 as the self-gravity operators which are 
defined below: 
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The coefficients of the operator O g are given below: 
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And we call Op and O 2 as the pressure operators which are 
defined below: 



The coefficients of the operator Op are given below: 
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The other pressure operator Op is given by: 
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To proceed further with the above calculations we need to 
know <jr and a v for the disc. We use the axisymmetric lo- 
cal stability parameter Q (Toomre, 1964) to get the radial 
velocity dispersion in the disc: 

ctrk 



Q = 



(16) 



3.36GE 

And for the azimuthal velocity dispersion we use the relation 
(Jip/cTR = k/2Q (Binney & Tremaine, 1987). For simplicity 
we consider a constant Q value for the disc. Note that by 
considering a constant Q we are making the perturbed pres- 
sures arbitrary by that constant. Since we do not have an 
apriori knowledge about the components of the stress tensor 
(hence a correct pressure) from the equations of motion, we 
are justified in playing around with the constant Q to eva- 
lute the perturbed pressures for our analysis. Our analysis 
for a cold disc (Q=0) remains unaffected by such assump- 
tion. A proper variation of Q on the radial co-ordinate will 
be considered in a later work. 

Now we require the perturbed surface density (E') and the 
potential ($') to be connected through the Poisson equation 
in order to produce a self-consistent solution of the prob- 
lem (eq.[8]). This we achieve using the integral form of the 
Poisson equation for the perturbed disc under the imposed 
m — 1 lopsided mode: 

$'(i?) = -G dR'R'Hi p(R,R')T.'(R') (17) 



Where the kernel in eq.(17) is given by 



Hi op (R, R') 



f 

Jo 



[R 2 + R' 2 - 2RR' cos a + 6 2 ] 2 



cos ada R 
n _ 

(18) 

The kernel represents the softened self-gravity of the per- 
turbation, 6 being the softening parameter. This allows us 
to perform the numerical integration over the nearby rings 
by removing the singularity at R = R' . We will discuss the 
effect of this softening parameter as we progress. The sec- 
ond indirect term (Papaloizou, 2002) in the above kernel 
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arises due to the m = 1 lopsided mass distribution about 
the geometrical centre of the disc. The indirect term plays 
a crucial role in making the disc susceptible to the lopsided 
instability. The importance of the indirect term is discussed 
in various places in the paper. 

Substituting the perturbed potential (eq.[17]) into 
eq.(8) we can arive at the following equation: 

cj 2 E'+£> iop (£> + S iop (£') = (19) 

where 

roc 

S iop (E') = -2lu p Z' + G / dR'R'Ci(R,R')T,'(R') 
Jo 

-OpE' (20a) 



Si op (E') =WpE' -u p G / dR'R'Ci(R, P')E'(P') 
Jo 
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+w p c£e' + eG / dR'R'C 2 (R, R')X'(R') - e0 2 p T,' (20b) 
Jo 

The kernels in the integrals of the above eq.(20a) 
and eq.(20b) are given by d(R,R') = 0\Hi p(R, R') and 
C2(R,R') = OgH.i op (R, R'). 

The above integro-differential equation (eq.[19]) describes 
the global behaviour of the m = 1 lopsided mode un- 
der the collective effects of self-gravity and pressure in the 
disc. Eq.(19) can be solved by recasting it into a matrix- 
eigenvalue problem. By discretizing on a uniform grid with 
N radial points in the disc we can write eq.(19) in a compact 
form: 

[lufop + DiopUJiop + Slop] Ei p = (21) 

Where I, T>i op , and Si op are the three N x N real square 
matrices. We call ®i op as the general damping matrix and 
Si p as the stiffness matrix for the galactic disc divided into 
N concentric rings in analogy with the nomenclature of lin- 
ear mechanical systems. Note that the damping, in our dis- 
crete N ring system i.e. in the disc, can either be positive 
or negative in nature. Negative damping will lead to a self- 
excitation of the discrete normal modes in the disc. In the 
above equation E; op is the eigenvector corresponding to the 
eigenvalue ui 

op • 

The matrix elements are evaluated below: 

Iij = Sij 

®Z P = -ZSiiMRj) + Mi(Ri, Rj) - Pi (ft, Rj) (22a) 



Slip = kjullRj) - u; p (R i )Mi(Ri,R j )+cj p (Ri)Pi(Ri,Rj) 

+e(R t )M2(Ri,Rj) - e(R i )P 2 (Ri,R j ) (22b) 

Where nMi(Ri,Rj) = ARRjCi (R t , Rj) and 
2TYM 2 (Ri,Rj) = ARR j C 2 (R l ,R j ) We use the method 
of finite difference to evaluate the contribution to the 
matrix elements from the pressure terms e.g. Pi (Ri, Rj) and 
P 2 (R l ,R j ). 



2.1 Technique to solve the equations 

Eq.(21) represents an TV- dimensional quadratic eigenvalue 
problem(QEP) in the eigenvalue u)i op and JV-dimensional 
eigenvector E ;o p describing the lopsided perturbation sur- 
face density in the disc. We consider N = 100 rings uni- 
formly spaced up to an outer boundary R ou t which we keep 
as a free parameter in our study. Since the original eq.(19) is 
an integro-differential equation we need to supply boundary 
conditions to solve the matrix problem. We use the Neu- 
mann boundary conditions dE'/dP = on the boundary in 
order to evaluate the contributions to the matrix elements 
from the pressure terms. Note that in the absence of veloc- 
ity dispersion i.e. for a cold disc we have a simple integral 
equations. The above boundary condition is used only to 
solve the govering equation in a hot disc. For a cold disc, we 
do not need to impose any external boundary condition; the 
boundary condition is in-built in the integral equation. The 
standard way to solve the QEP (eq. [21]) is to reduce it to 
a generalized eigenvalue problem(GEP). We use LAPACK 
subroutines for the numerical solution of our GEP. For more 
numerical details see Saha & Jog (2006) and relevant refer- 
ences therein. 



3 RESULTS 

We study the lopsided mode in a galactic exponential disc 
with the central surface density Eo and a scale- length Rd\ 

E°(P) = T, e- R/Rd (23) 

as observed in most spiral galaxies (Freeman 1970). We limit 
our investigations up to 4 disc scale- lengths, typically within 
the optical region, which contains ~ 90% of the total stellar 
disc mass. We do not include the effect of the dark matter 
halo in the present problem since different dynamical mech- 
anisms will then be involved, and we want to examine them 
in turn. 

In all the calculations we get the perturbation surface den- 
sity i.e. eigenvector Ei op (eq.[21]) in a dimensionless form 
and the eigenvalues u>i op corresponding to the eigenvector 
are in units of \J nGEo/Rd- For disc like our Milky Way the 
value of this unit ~ 52.9 kms _1 kpc _1 . 

We first consider here a cold self-gravitating disc i.e a disc 
with zero velocity dispersion. The unperturbed disc is rota- 
tionally supported against the gravitational attraction. The 
underlying motivation is to find out how the behaviour of 
lopsidedness (m = 1 mode) in a differentially rotating cold 
self-gravitating disc changes when a finite non-zero velocity 
dispersion is included in the disc. Below we report the re- 
sults first for a cold disc and subsequently the results from 
the hot disc are discussed. 



3.1 Lopsidedness in a Cold Disc : 

We show that the cold self-gravitating discs are able to sup- 
port m — 1 discrete normal lopsided modes. These modes 
are important because they show the observational signa- 
tures and they also turn out to be long-lived as discussed 
below. 
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Figure 1. Behaviour of lopsided modes in a cold disc. Mode 
shapes arc quite robust with respect to the changes in the disc 
size. The amplitudes of these modes arc in arbitrary units. 



3.1.1 Lopsided mass distribution in stellar disc 

In Fig. 1 we have shown four subplots of the lopsided 
modes at various values of the outer boundary Rout- The N- 
dimensional matrix eigenvalue equation (eq.[21]) produces 
2N eigenvalues and 2N eigenvectors in the disc. Observa- 
tionally relevant modes are actually the lowest frequency 
ground state modes. By ground state we mean the mode 
with the lowest value of the real part of the eigenvalue 
(Re(cj iop )) in the eigen spectrum. These modes are with the 
lowest number of nodes, typically zero or one node at most, 
showing the global nature of the lopsidedness in the disc. 
For smaller sized disc (as shown in the first two panels in 
fig. (1)) these modes are very interesting in nature, they are 
discrete and distinct (as can be seen in fig. (4)), in the sense 
that these modes represent a very slowly precessing station- 
ary lopsided pattern. The imaginary parts of the eigenvalues 
corresponding to these pattern are zero i.e. Im(o)j op )=0 for 
these modes. However as we increase the disc size this sce- 
nario changes and collective effect due to disc self-gravity 
becomes very important. We find that for a comparatively 
larger disc the most slowly precessing pattern appears as an 
instability in the system. The reason for this instability lies 
in the collective effect due to the softened self-gravity. Even 
though its hard to exactly pinpoint the source of the lop- 
sided instability, our numerical analysis show that without 
the indirect term the ground-state lopsided modes are never 
unstable in the disc. 

We have used the value of softening parameter b ~ some 
fraction of the inter ring spacings in order to solve the ma- 
trix eigenvalue problem. We show that the mode shape cor- 
responding to the lowest pattern speed is quite robust with 
respect to the variations in both the parameters namely R ou t 
and b. In Figs. la - Id, we have varied R ou t from 2.8Rd to 
3.8Rd and clearly in all these cases the mode amplitudes in- 
crease with the galactocentric radius, in a good agreement 
with the observations (Section 1). 

The surface density contours for the lopsided disc : 
T,(R,<fi,t) = T,°(R) + S p T,'(R,if,t) are shown in Fig. 2. We 
have considered 8 P = 0.02 and the contours are at t=0. This 
clearly brings out that the contours deviate from the unper- 
turbed disc surface density more at larger radii. A barycen- 
tric shift in the mass distribution of the disc due to the lop- 
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Figure 2. Surface density contours of the lopsided disc for S p = 
0.02. The maximum of the surface density occurs at (0,0). The x 
and y axes are in units of R d . The outward contours are more and 
more deviated from the unperturbed circular contours — showing 
that lopsided behaviour is prominent in the outward direction. 
The contour levels are 0.05, 0.1, ...0.3 times the central surface 
density. 



sided mode is obvious from fig. 2. Such a displacement of the 
center of mass of the disc due to a one-armed spiral mode is 
discussed in great detail in Evans & Reed, 1998. They show 
that a one-armed spiral does not shift the barycentre of the 
disc substantially. The lopsided mode in our analysis seems 
to precess with a very slow pattern speed (as can be seen 
in fig. 3 and fig. 8). Because of this slow pattern of oscillation 
we expect the barycentre of the disc to be shifted only by a 
small amount and the mass distribution of the disc does not 
alter appreciably; it probably stays in that state for a long 
time. However, since in the linear regime the amplitude of 
the lopsided surface density is arbitrary, we do not think it 
is possible to calculate the exact shift of the barycentre due 
to the lopsided pattern. 

So our study - based on the internal disc dynamics- nat- 
urally explains the observational fact that the lopsidedness 
increases outward in the disc. 



3.1.2 Persistence of the lopsided mode 

In order to achieve a persistent lopsided pattern in the disc, 
the pattern formed initially by the stellar orbits should be 
synchronized in such a way that they rotate with a con- 
stant pattern speed. In the kinematical picture this requires 
that the differential precession in Q — k, must be zero which 
is true in a Keplerian disc because the disc self-gravity is 
insignificant and the lopsided pattern would remain intact 
there for all time. In a sharp contrast, in a galactic stellar 
disc, it is the self-gravity which is the most important. So 
the survival of lopsidedness in the stellar disc is puzzling. 
The strong differential precession in n — f2 in the inner re- 
gion would wash away any coherent pattern imposed onto it. 
The typical winding up time scale for the lopsided pattern 
can be obtained by using the formula r iop = 2-k/A(k — fl) 
(See BLS). Their estimation shows that the typical lifetime 
of the lopsided pattern seen in outer HI disc is ~ 1 — 5 Gyr. 
For a typical stellar disc like that of our Galaxy this winding 
up time scale in the optical region would be less than a Gyr. 
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Figure 3. The solid line, connecting the open circles, shows the 
precession rate Re(uii op ) of the lopsided mode under the collective 
effect of the self-gravity. The dashed line is free precession and 
acting as the lower boundary (Q — k) of the continuum. The upper 
boundary (Q + k) is not shown because it is orders of magnitude 
higher than Re(u>i op ). 



So in the stellar disc the survival problem of the lopsidedness 
is far more severe. 

In Fig. 3 we have shown the variation of the Re{tui op ), 
denoting the pattern speed of the lopsided mode in the expo- 
nential disc, with respect to the size of the disc. The eigenfre- 
quencies uii op are corresponding to the ground-state lopsided 
mode in the disc. We see that Re{uii op ) is almost constant 
compared to the variation in 57 — k and thereby the disc 
avoids the differential precession. In this way the winding- 
up time scale increases by a factor of ~ 10 or more and this 
turns out to be more than 12 Gyr for a disc like that of our 
Galaxy. This implies that the lowest frequency modes are 
almost non-rotating pattern of the disc and N-body simula- 
tions (as described later in sec. 4) confirms this picture. Note 
that in all our calculations Re(u)i op ) lie in between the gap 
made by Q — k and Q. + k which act as the lower and up- 
per boundary of the continuum region respectively. In other 
words we have 57 — k < Re(uoi op ) < Q + k. All the ground- 
state lopsided modes in our calculations are thus discrete in 
nature; once Re(uji op ) lies outside this gap they fall in the 
continuum regime. We show that as the disc size increases 
the width of the gap shrinks and eventually any discrete lop- 
sided mode ceases to exist! And this phenomenon occurs at 
a region near 4.6 disc scale-lengths where Re(iOi op ) ~ n — Q 
(see fig. 3). Surprisingly there is a natural resonance at k — Q 
at around 4.6 disc scale-lengths. So we see that vanishing of 
the discrete lopsided modes occurs below this resonance and 
this might be acting naturally as a boundary beyond which 
no discrete lopsided mode can exist in a purely exponential 
disc. 

Fig. 3 also brings out another fact clearly regarding the sense 
of precession of the lopsided pattern in the cold disc. In lin- 
ear analysis the final equation (21) is invarinat under the 
transformation uJi op — > — cD iop , one would expect the pattern 
speed to be both retrograde and prograde. Our numerical 
eigenvalue analysis shows that both the retrograde and pro- 
grade pattern are likely to occur in the cold disc, some of the 
earlier works namely Statler (1999) has argued that the pre- 
cession of a self-consistent lopsided disc must be prograde; 



Jacobs & Sellwood (2001) find in their simulation that the 
lopsided modes precess in prograde direction. 



3.1.3 Complete eigen spectrum of the lopsided 
pattern in a cold disc 

On solving the N-dimensional eigenvalue problem (eq.[21]) 
in the case of cold discs (Q=0) with different disc sizes, we 
obtain 2N complex eigenvalues and these complex eigenval- 
ues are plotted in the Argand diagram as shown in fig. 4. The 
complete eigen spectrum in the case of spiral and bar-like 
structures in galaxies are discussed by Polyachenko, 2004. 
The complex eigenvalues appear as complex conjugate pairs 
in the full spectrum. In the Argand diagram, we will concen- 
trate only to those points for which Re(o;; p) ~ since we 
are primarily interested in the slow modes. For smaller sized 
disc as in panel-a, there exists a distinct (clearly separated 
from the rest of the points in the complex plane) point in the 
discrete region close to the boundary of the lower continuum 
(denoted by the vertical solid line) corresponding to Re(w; op ) 
~ and lm(cui op )—0. The mode corresponding to this point 
is an almost non-rotating stationary global lopsided pattern 
of the disc. Being isolated in the Argand diagram this mode 
resembles like an excited normal mode of the disc much like 
what Lynden-Bell in 1965 first conjectured in the context of 
warps (basically m = 1 mode in the perpendicular direction 
to the disc midplane) that they are excited normal modes of 
the discs. The existence of such a mode of oscillation outside 
the continuum was shown by Mathur, 1990 for one dimen- 
sional and 3D gravitating system. The presence of such a 
distinct normal mode in our eigen solution for a disc-like 
2D self- gravitating system confirms the previous trend found 
in theoretical studies namely self-gravitating system supports 
normal modes of oscillation in the disc. 

Interestingly, note the appearance of a complex conju- 
gate pair of eigenvalues in the continuum region near the 
solid line in panel-a. As the disc size increases this pair 
drifts towards the discrete region from the continuum and at 
around 4 disc scale lengths, they are in the discrete region. 
And the slowly precessing stationary pattern turns out to 
be a lopsided instablity in the cold disc. This instability, in 
comparatively larger discs, arises due to the collective effect 
of the softened self-gravity of the perturbation where the 
indirect term place a crucial role in driving the instability. 
Whereas for smaller sized discs, which are like a broad an- 
nular ring, the collective effect of softened self-gravity is not 
enough to drive an instability and the lopsided pattern ba- 
sically represents a stationary oscillation of the whole disc. 
Apart from the discrete region, the complete eigen spectrum 
contains a lot of modes in the continuum region. For a cold 
disc only a few of these modes are with lm(a; ;op )=0 repre- 
senting stationary oscillation in the system. These particu- 
lar normal modes (as shown in fig. 5) are interesting because 
they resemble the van Kampen modes in plasma physics. 
These modes would probably decay by Landau damping and 
hence are of no interest in our present study. 



3.2 Lopsidedness in a hot disc : 

In reality galactic discs are never cold. Because of a non- 
zero velocity dispersion a hot disc is supported by rotation 
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Figure 4. A plot showing the behaviour of the complex eigen- 
values corresponding to the lopsided modes in a cold disc (Q=0). 
Each panel shows the complete cigen solution of cq. [21] for a dif- 
ferent value of the disc size R ou t- The vertical dashed line is the 
lower boundary of the continuum as shown in fig. 3 for each disc 
size. 
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Figure 5. The plot showing the typical behaviour of a stationary 
van Kampen mode lying in the continuum for Rout=3.0. These 
continuum modes are corresponding to lm(aj; of> )=0. 



as well as pressure against the gravitation. We quantify the 
hotness of the discs by the Toomre's Q (eq.[16]) parame- 
ter; note that Q=0 refers to a cold disc and Q/0 means a 
hot disc. If Q is less than the critical value of 1, the disc is 
axisymmetrically unstable. Basically Q serves as a very use- 
ful local criterion for the axisymmetric stability of the disc. 
Normally Q has a radial profile in the disc. For the sake of 
simplicity we consider Q=constant throughout the disc. By 
doing this we would like to understand what happens to the 
lopsided perturbations when there is a non-zero Q in the 
disc. 



3.2.1 Complete eigen spectrum of the lopsided 
pattern in a hot disc 

The behaviour of the cold disc against lopsided perturbation 
changes dramatically as soon as a finite velocity dispersion 
is introduced in the system. The transformation in the disc 
behaviour is best seen through the complete eigen spectrum 
plotted (fig.6) in the Argand diagram for different values of 



Q. In the case of Q=0, there are only a few modes in the con- 
tinuum with Im(ci;; O p)=0 and only one or two in the discrete 
region. Most of the modes in the continuum are of unsta- 
ble in nature. However as we increase the Q value, the disc 
changes abruptly and modes in the continuum are populated 
along the real line i.e. with Im(a;; O p)=0. For Q=0.2 in this 
case shows that there are almost no unstable modes in the 
discrete region, some of them are indeed in the continuum 
with high pattern speeds. At Q=0.6, we see the appearence 
of unstable modes in the discrete region of the spectrum. 
The appearence of the unstable mode in the spectrum is not 
very obvious. There are subtle interplay amongst the indi- 
rect term in the poisson equation (eq.[18]), the softened part 
of the self-gravity and the velocity dispersion (Q). Without 
the indirect term the ground state modes are never unsta- 
ble, no matter what the values of the softening (b) and the 
Q-value in the disc. The indirect term turns out to be the 
main driver of the lopsided instability in the disc. 

Beyond this as Q increases, the number of unstable 
modes goes down and at some Q=Q C rit no unstable modes 
are found with lowest pattern speed. Beyond Qcrit the 
modes with lowest value of the pattern speed in hot disc 
again turn out to be purely oscillatory. 



3.2.2 Dependence of the growth rate on velocity 
dispersion 

We solve the governing eq.(21) for the lopsided modes in a 
hot disc with non-zero radial (ctr) and azimuthal (u v ) veloc- 
ity dispersions. Inclusion of velocity dispersion has an effect 
opposite to that of the self-gravity in the disc. In a cold disc, 
the imposed lopsided perturbation emerges to be a station- 
ary normal mode oscillation for smaller sized disc or grows 
as an instability when the collective effect of the softened 
self-gravity is sufficient as it is in comparatively larger sized 
discs. By addition of a finite non-zero velocity dispersion 
(quantified here by a constant value of the Q-parameter) in 
the disc the strong effect of self-gravity is diluted. As we 
increase the Q-value the growth rate (7; p=Im((x>; p)) de- 
creases and there exists a Q cr ; t beyond which we find that 
the lowest frequency modes are no longer unstable in the 
disc, see Fig. 7. Note that as we increase the Q-value, the 
growth rate reaches a maximum and starts to decline; a 
similar behaviour of the growth rate vs Q was observed in 
the global analysis of m = 1 mode in a nearly Keplerian disc 
(Adams et al. 1989). In Fig. 7 we find the value of Q cr it to 
be 1.0 for a disc with R ou t=3.5 R<j, for other values of Rout 
we find that Q cr it normally lies around the 1. The value 
of Qcrit may be surprising to the reader; note that we are 
considering here a purely stellar disc without any dark mat- 
ter halo. Take for example the stellar disc of our Galaxy 
with central surface density So = 640 M0pc~ 2 , 11^=3.2 kpc 
and the radial velocity dispersion (<tr) as according to the 
Lewis-Freeman law (Lewis & Freeman 1989), the Toomre's 
Q parameter is always less than 1 beyond 1 disc scale length. 



3.2.3 Dependence of the lopsided pattern on 
velocity dispersion 

In Fig. 9 we study the behavior of the lopsided modes in a hot 
disc. These lopsided modes are again the lowest frequency 
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Figure 6. The full eigen spectrum for a disc at different val- 
ues of Q, including one with Q=0.0. The figure shows how the 
eigen spectrum behaves in the Argand diagram as finite velocity 
dispersion is introduced in the cold disc of size Rout = 3.5-Rd- 
The verticval dashed line represents the lower boundary of the 
continuum as shown in fig. 3 for the disc. 
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Figure 7. Growth rates (7; op ) of the lowest frequency modes 
are plotted against the Q-parameter. The growth rate reduces as 
Q increases in the disc followed by a peak. Beyond Q cr ;t these 
lopsided modes are not unstable anymore. 



(i.e. the lowest value of the real part of the complex eigen- 
value) eigenmodes out of the 2N eigenmodes in the N-ring 
eigen-system. In a cold disc these lowest frequency eigen- 
modes are interesting or have observational relevance be- 
cause they are global in extent (having zero or one radial 
node). With the introduction of a finite velocity dispersion, 
the behaviour of the eigenmodes also changes drastically. 
In Fig. 9 the solid line is the eigen mode in a cold disc i.e. 
Q=0. As Q-value increases in the disc, the self-gravity starts 
falling down and fails to synchronize the adjacent rings to 
precess them coherently and produce a global mode. The 
eigen modes start acquiring radial nodes. The number of 
radial nodes increases as Q approaches the value Q cr it- At 
around Qcrit the eigen modes have vanishing amplitudes 
with many radial nodes beyond 1 disc scale length show- 
ing that there is no coherent lopsided pattern anymore in 
the disc beyond Q=Q cr it- In fig. 8 we show the variation of 
the pattern speed with the Q-value. Interestingly the pat- 
tern speed does not get affected in appreciable amount with 
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Figure 8. A plot of the pattern speed vs the Q value in the disc 
for a disc size of 3.5 Rd- Clearly, there is no preferred sense of 
precession of the lopsided pattern in the disc. 
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Figure 9. The lowest frequency modes are plotted against the 
radius for different values of the Q-parameters at t=0. As Q in- 
creases in the disc, the number of radial nodes of these modes 
increases. The solid curve is for Q=0.2; dashed one is for Q=0.6 
and the dotted one is for Q cr it=1.0. 

. Note that when Q=0.0 i.e. for a cold disc, the lopside modes are 
coherent with zero or one node at most ( fig. 1). This coherence 
in the lopsidedness is lost as the disc becomes hotter. 

the increase in Q. The pattern speeds of the ground state 
modes are not quite insensitive to the increment in the ve- 
locity dispersion. Our numerical calculation show that there 
is no definite sense of precession of the lopsided pattern in 
the disc. Fig. 8 also brings out clearly another fact that the 
global lopsided pattern precesses quite slowly even in a hot 
disc! Even the dynamical friction (Ideta, 2002) from an ex- 
ternal dark matter halo would take long time to damp out 
such a pattern in the galactic disc. 



4 N-BODY SIMULATIONS 

In order to test the predictions of the previous analysis, 
numerical simulations have been performed. Since differ- 
ent physical mechanisms can occur in the presence of dark 
matter halo, and the presence of a gaseous component, the 



present work is restricted to a pure stellar disc. Initially there 
is no bulge component, although a pseudo- bulge is building 
along the simulations, through bar instability, which then 
triggers a box-peanut shaped bulge formation. 

These simulations are meant to isolate the dynamical 
mechanisms occurring in a pure stellar disc, for academic 
purposes. In a following paper, more realistic galaxies with 
dark matter haloes, and gaseous components will be consid- 
ered. 

4.1 Numerical details 

The galaxy is isolated, and represented by only stellar par- 
ticles. The most efficient N-body algorithm to simulate such 
a distribution is a PM (particle-mesh) code. Since in par- 
ticular m — 1 modes and lopsidedness are of interest, polar 
grids with higher spatial resolution in the center are not 
adapted: the high density core of the galaxy is expected to 
move around, and the variable resolution could produce arte- 
facts. A cartesian grid is therefore selected. We use here an 
FFT N-body code with the James (1977) method to avoid 
the influence of Fourier images. The useful grid is 256 2 xl28, 
corresponding to (60 kpc) 2 x30 kpc, implying a cell size of 
234 pc. This size is also the softening length of the Newto- 
nian gravity. The total number of particles is between 8 10 s 
and 1.6 10 6 , according to the various runs. 

The initial radial distribution of particles is not very im- 
portant, since it is known that a pure disc without any spher- 
ical component (such as bulge or halo) is unstable with re- 
spect to bar formation; then the bar re-distributes efficiently 
the particles, to form an exponential disc (e.g. Combes et al 
1990, Pfenniger & Friedli 1991). For the sake of simplicity, 
the initial distribution of stellar particles is taken according 
to a Kuzmin-Toomre disc, of surface density: 

E(r) = E (l + r 2 /ri)- 3/2 

where is the characteristic scale-size of the disc, truncated 
at 25kpc, with a disc mass M^. The initial velocity distribu- 
tion of the stars in the disc is computed from the stationary 
Jeans equation in cylindrical coordinates, yielding the asym- 
metric drift Vrot — v C i r , assuming that a v /o r — k/(2Q), as 
in the epicyclic approximation. 

The distribution of density perpendicular to the disc 
is that of an isothermal thin layer: p(z) — posech 2 (zj H) 
with a height slightly flaring with radius, as H(r) = Hq(1 + 
(r/9fcpc) 2 ) 0,35 , where Ho = 700pc. The z-velocity dispersion 
is then computed from a 2 — H(r)irGYl(r). 

The initial Toomre Q parameter is slightly varying with 
radius as: Q = a r /a crit = Qo/exp[(rd/23kpc) 2 } The time 
step is 1 Myr. The parameters of the runs described here 
are displayed in Table [T] 

The stellar rotation curve corresponding to the refer- 
ence run A is plotted at the end of the simulation in Fig. 1101 
The characteristic frequencies fl, Q, — k/2 and SI + k/2 are 
derived from the measured potential, at the final epoch. 

4.2 Results 

The behaviour of run A and run B is quite similar. Run B 
has twice more particles, and demonstrates how the particle 
noise may affect the estimation of the m = 1 amplitudes 
(error bars will be given below). 
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Table 1. Parameters of the initial conditions 
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Figure 10. Rotation velocity Vrot of the stars, at the end of 
simulations (run A, T=12.6 Gyr), and frequencies H, Q — k/2 
and Q + k/2 derived from the potential at this epoch for run A. 

Our reference run (A or B) reveals a bar instability 
as expected, given that no dark matter halo, nor bulge is 
present to stabilise the exponential disc. The strength of the 
bar is estimated by the ratio of the tangential force of the 
m — 2 Fourier component, normalized to the radial force. If 
the potential is decomposed as 

$(r, tp) = $o(r) + ^2 $ m (r) cos(rrup - <j> m ) (1) 

m 

the bar strength is the maximum over radius of 

Qa = 2$ 3 /r|F r | (2) 

where F r — —d<&a/dr. The evolution of Q2 for run B is 
shown in Fig II 11 as well as its pattern speed fib, measured 
from the monitoring of the phase of the m = 2 pattern with 
time. 

In the beginning of the evolution, the m — 2 strength 
has a contribution from the spiral structure which develops, 
and helps to transfer angular momentum in the outer parts. 
After a few Gyrs, and in the quasi-stationary state, there is 
only a bar. 

The bar instability heats the stellar disc, which in a 
first phase weakens the m — 2 pattern that developped its 
maximum around lGyr. Then the bar continues to grow 
slowly until its strength reaches a stationary value, of Q2 ~ 
0.1. The bar pattern speed increases slowly in this evolution, 
which is due to the strong concentration of matter driven by 
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Figure 11. Strength of the m = 2 mode (Q2) and pattern speed 
£1;,, for run B. Units of Of, are lOOkm/s/kpc. Also shown is the 
m = 1 strength Q\. 



Figure 13. Radial behaviour of A\ for all runs, at the epoch 
T= 7.2 Gyr. A\ is the m = 1 Fourier component of the surface 
density, normalised with the m = component. 



the bar. The face-on morphology of the galaxy is displayed 
in the contours of Fig [12] The edge-on projections reveal 
the build-up of a peanut morphology, which considerably 
thickens the stellar disc, and forms a pseudo-bulge. 

The abrupt fall-off of the Q2 strength in the first evolu- 
tion phases (after 2-3 Gyr), is due to the violent instability 
of the beginning, occuring so shortly with respect to the dy- 
namical time, that the family of orbits sustaining the bar 
have no time to re-arrange, before the radial redistribution 
of matter. The mass is re-distributed radially, in consequence 
of the angular momentum flow towards the outskirts. The 
exponential disk length decreases, and the concentration of 
matter has the consequence that the bar pattern speed in- 
creases. Then the resonances change radial location, and 
therefore involve different particles, which means that the 
bar can now grow again with a different pattern speed, and 
different sustaining orbits. 

To study the simultaneous development of any m = 
1 mode, the analogous Qi term was computed from the 
Fourier analysis of the potential. However, the m = 1 ampli- 
tude is weaker, and the estimator from the potential, which 
is a more global way to quantify the lopsidedness, dilutes any 
local maximum, and therefore revealed too noisy. An estima- 
tion from the density map in the face-on projection of the 
galaxy, was then preferred. This is an estimator currently 
chosen by observers (e.g. Zaritsky & Rix 1997, Bournaud et 
al 2005). From the Fourier analysis of the surface density in 
the plane, we define Ai the ratio of the m = 1 and m = 
term. The behaviour of A\ as a function of radius is dis- 
played for all runs in Fig 1131 For run C and D, there is a 
strong maximum in the center of the disk (around 5 kpc). 
After a pronounced minimum around 14kpc, the amplitude 
of Al is rising again in the outer parts. The projected density 
maps of the perturbed disks show that this m = 1 is essen- 
tially due to the off-centering of the corresponding parts of 
the disk. The central parts are shifted on one side, while the 
outers parts are shifted on the opposite side. This explains 
the minimum in the intermediate radii (around 14kpc). Al- 
though there is mainly at the beginning an m = 1 spiral 



structure, the dominating feature of the m — 1 mode is the 
off-centering. 

It can be noted that the estimator gives a coherent value 
getting out of the noise between 2 and lOkpc. Values at large 
radii are always noisier, due to the low number of particles, 
in the outer parts of an exponentially decreasing surface 
density. In the very center, the small area involved also does 
not include enough particles. We therefore select for all runs 
the common radius R= 5kpc to estimate Ai and deduce the 
time variations. 

While Run B is unstable to the m — 2 mode (bar and 
spiral), it does not show a strong m = 1 perturbation, which 
could be expected, given the high initial velocity dispersion 
(Q = 2.3). Run C on the contrary is violently unstable with 
respect to m = 1, as can be seen in the evolution of pattern 
strength in Fig [15] and in the contour plots of Figure 1161 
This is easily explained through the lower initial velocity 
dispersion Q = 1.2. 

Both runs develop a bar however, and are unstable with 
respect to the buckling instability, forming the now well- 
known peanut shape (e.g. Combes et al 1990, Martinez- 
Valpuesta et al 2006). The peanut formation interferes with 
the m = 1 mode to make the z-elevation of particles asym- 
metric, and produce a marked tilt in the disc plane at the end 
of the simulations. This does not occur for run B, where the 
distribution of particles remain roughly symmetrical with 
respect to the center of mass. But Run C shows very well 
this phenomenon. In the buckling some particles are ele- 
vated upwards, at some radii, in resonance with the vertical 
Lindblad resonance, which happens to coincide with the in- 
plane ILR, since k ~ u z , at these radii. The lopsidedness 
happens to give more weight to the population of particles 
that are elevated upward, breaking the symmetry with re- 
spect to the plane. Some particles escape at larger radii, and 
are lost for the equilibrium. The remaining disc re-settle in 
a tilted orientation. 

Bending instabilities in galactic discs are expected when 
the ratio between z and r velocity dispersion (cr z /oy) is lower 
than the Araki (1985)'s value of 0.3 (e.g. Revaz & Pfenniger 
2004). In our runs, this ratio is initially about 0.5 in the 
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Figure 12. Run B: Contours in logarithmic scale of the surface density of the stellar disc, face-on (top) and edge-on (bottom), at 
four different epochs: T=0.8, 5.6, 9.6, 14.4 Gyr, from left to right. The contours are taken between the maximum surface density, and a 
minimum equal to 5 10 — 4 that value. In units of 900 Mq/pc 2 , in square pixels of 470pc in size, the maxima are respectively 0.6,0.6,0.6,0.6 
for the face-on values, and 1.6, 1.9, 1.5 and 1.5 for the edge-on values. 




center, and falls towards 0.1 in the outer disc. Instabilities 
are therefore expected. At the end of the simulation, the 
a- z /oy ratio increases to be just above 0.3 (cf Figure [T4)) . 

In the edge-on views of Fig. 1161 it is easy to see the 
shape of the tilt (see the T= 9.6 Gyr epoch for instance). 
The center of the disc until a radius of llkpc is tilted with 
a position angle larger than 90°, while the outer parts are 
tilted the opposite way (with position angle smaller than 
90°). The disc shows a warp, with a line of nodes following 
the peak of the m = 1 perturbation. The line of nodes is 
almost straight. The Ai amplitude is also conspicuous inside 
and outside a radius of about 14kpc, as is visible in Fig. 
1131 This is due to the off-centring of the central part which 
is opposite of that of the outer parts. The limiting radius 



corresponds roughly to the resonant radius where £1 — k — 
0. 

Run D was carried out with a twice smaller exponen- 
tial scale-length, in order to explore the varying smoothing 
length parameter. The latter is equal to the cell size of 234 
pc, so that the ratio of the disc characteristic size to the 
softening length is reduced by a factor 2. The resulting disc 
shows a larger instability than for the reference run B, as 
can can be seen in the evolution of pattern strength in Fig 
[T71 and in the contour plots of Figure [TS] 

At first, it may appear paradoxical that run D, charac- 
terized by a larger softening with respect to the disc scale, 
is more unstable to the global m = 1 mode. For usual WKB 
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Figure 16. Run C: Contours in logarithmic scale of the surface density of the stellar disc, face-on (top) and edge-on (bottom), at 
four different epochs: T=0.8, 5.6, 9.6, 14.4 Gyr, from left to right. The contours are taken between the maximum surface density, and 
a minimum equal to 5 10 -4 that value. In units of 900 Mq/pc 2 , in square pixels of 470pc in size, the maxima are respectively 1.2, 1.7, 
1.8, 1.9 for the face-on values, and 3.6, 2.5, 2.5 and 3.2 for the edge-on values. 
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Figure 17. Bar strength (Q2) and bar pattern speed fit,, for 
run D. Units of £7^ are lOOkm/s/kpc. Also shown is the m = 1 
strength Q±. 



Figure 19. Strength of the m = 1 mode, for run B, C and D. A\ 
is the m = 1 Fourier component of the surface density, normalised 
with the m = component, and taken at R=5kpc. 



local waves, self-gravity is reduced by the softening, and disc 
should be more stable. But for these global modes, involv- 
ing the whole disc, with a wavelength comparable to the disc 
size, this is different. The instabilities are due only to the 
collective effect of gravity, and the softening mimicks the ef- 
fect of pressure, which stabilises more local modes. The disc 
is then colder, and can sustain global modes. 

The time evolution of the m = 1 mode is compared for 
all runs in Fig |19l where is plotted the value of A\ estimated 
at the radius of 5kpc. A maximum at a value of A\ ~ 2 
is reached around T= 7 Gyr for run C, and at about the 
same time for run D, at A\ ~ 1. For run B the value is 
negligible and tends to be stationary at ~ 0.1. This evolution 



corresponds to a coherent lopsidedness deformation of the 
stellar disc. The deformation appears quasi still in phase, 
and within the noise it is impossible to distinguish a value 
for the pattern speed fh different from zero. 

The m = 1 perturbations found in the present sim- 
ulations might appear quite large, with respect to what is 
usually found in galactic simulations. However, this is due to 
the fact that in general galaxy models are assumed to be em- 
bedded in massive spherical haloes, which stabilise the discs. 
Here we have purely stellar discs, initially without any cen- 
tral bulge component. Simulations with massive discs have 
also been run by Revaz & Pfenniger (2004), who also find 
m = 1 perturbations quite common. 
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Figure 18. Run D: Contours in logarithmic scale of the surface density of the stellar disc, face-on (top) and edge-on (bottom), at 
four different epochs: T=0.8, 5.6, 9.6, 14.4 Gyr, from left to right. The contours are taken between the maximum surface density, and 
a minimum equal to 5 10 -4 that value. In units of 900 Mq/pc 2 , in square pixels of 470pc in size, the maxima are respectively 1.7, 1.7, 
1.7, 1.7 for the face-on values, and 4.4, 3.8, 4.0 and 3.9 for the edge-on values. 



5 COMPARISON BETWEEN ANALYTICAL 
ANALYSIS AND NUMERICAL 
SIMULATIONS 

The numerical simulations of an isolated purely stellar ex- 
ponential disc reveals some confirmation of the analytical 
derivations about the global mode. Although the m — 2 
perturbation is always very strong, which is a feature not 
considered in the calculations, it is possible to see the devel- 
opment of an m = 1 perturbation, which appears to involve 
the whole disc, with a global characteristic. The wavelength 
is long, at least of the same size as the disc. 

The instability in the simulations appear for values of 
the Toomre parameter larger than is predicted for the ana- 
lytical calculations. But the latter are computed for an in- 
finitely thin disc, and do not take into account all physical 
parameters, except with approximations. The "effective" Q 
in the analytical model should be larger than computed, to 
take account of the pressure, the softening, the thickness of 
the plane, for instance. Merritt & Sellwood (1994) found also 
that large-scale modes could develop in numerical galaxies, 
well beyond the predictions of the analytical calculations in 
a thin sheet. 

All runs develop a bar, a strong m = 2 perturbation, 
that slowly develops a peanut instability. The disc is buck- 
ling and thickening around the inner Lindblad resonances, 
between 4 and 6 kpc in radius for runs B and C (3kpc for 
run D). The disc thickening is also enhanced by the heating 
of the disc through the gravitational instabilities, and sat- 
urates after 10 Gyr, as shown in Fig. 1201 The thickening of 
the plane is not only enhanced by the peanut alone (as in 
run B), but also by its coupling with the m = 1, leading to 
the plane tilt: the final thickness is higher in run C where 
both are present (Fig. I20|) . 

For a given softening length, the m = 1 develops es- 
sentially for the initially colder discs, i.e. initial low values 
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Figure 20. Average thickness of the plane (variance of the z- 
coordinate averaged over all particles) for run B, C and D. The 
time evolution reveals a saturation after 10 Gyr. 



of the Toomre parameter Q. The run B develops mainly a 
bar, with little asymmetries, practically no m = 1, while 
the initially cold disc of run C develops a pronounced lop- 
sidedness. This strong instability heats violently the disc, 
which is much hotter in the center for run C, as shown by 
the Toomre parameter averaged over the radii inside 5 kpc 
(Fig. I2ip . Since the disc is lopsided, the exact calculation 
of this parameter is not possible, the center of mass of the 
galaxy being different from the maximum of density, or the 
minimum of the gravitational potential. However, averaged 
over a much larger radius, the heating of the all runs become 
more similar (Fig. I22[) . 

The pattern speed of the m — 1 is very slow, such that 
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Figure 21. Evolution of the Toomre parameter, for all runs, 
averaged over the center of the disc, for radii lower than 5 kpc. 
The average is weighted by the number of particles at each radius. 



Figure 22. Evolution of the Toomre parameter, averaged for 
radii lower than 10 kpc. The saturation around Q = 2 is then 
visible for all runs. 



it is difficult to measure. The perturbation is almost still, 
difficult to estimate as prograde or retrograde, but certainly 
is partly retrograde at large radii. This confirms our semi- 
analytical results about the pattern speed of the lopsided 
mode. 

A more interesting feature is the behaviour with respect 
to the relative softening length b (ratio of the softening to 
the exponential disc scale). The instability with respect to 
the m — 1 global mode is larger for a larger b (run D), than 
for a smaller b (run B), at the same value of Q. This can be 
interpreted in terms of global mode due only to collective 
gravity effects. 

This behaviour might appear surprising, as the soften- 
ing is known to stabilize discs with respect to gravitational 
instabilities. However, the softening suppresses essentially 
the small scale perturbations, the modes with small wave- 
length (obeying the WKB approximation) . This then leaves 
a disc more receptive to global modes, and instabilities with 
wavelength of the order of the disc radius. 

More quantitatively, Romeo (1994) has shown that the 
softening s weakens the potential perturbations of radial 
wavenumber k — 2ir/\, by a factor exp(—\k\s). This effi- 
ciently stabilises all small scale instabilities for run D, which 
is then colder in the first Gyr of evolution and more able to 
develop global modes. 

Another argument that could explain that large wave- 
length modes are favored in runD with respect to run B, is 
the order of magnitude of the self-gravitating critical scale 
length, X crit = 4tt 2 G£/k 2 (cf Toomre 1981). Since the scale 
length rd of the disc in run D is divided by 2, keeping the 
total mass the same, the surface density is multiplied by 4, 
and k 2 is multiplied by 8. The critical scale \ cr it is divided 
by 2. However, all initial discs were truncated at R=25kpc, 
and the disc extent is larger with respect to the critical scale 
in run D. 

The efficiency of the swing amplification is parametrized 
by the X factor (Toomre 1981) where X — A/(A C r«sini), and 
i is the pitch angle, nearly 90° here. For wavelengths of the 
order of the disc extent, this X factor is initially larger for 
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Figure 23. Evolution of the X parameter, indicator of the swing 
amplification efficiency, averaged for radii between 10 and 20 kpc. 



run D, as shown in Fig 1231 and approaches the optimum 
value for maximum amplification, which is 2. 

As in the calculations, the simulations show that the 
lopsided perturbation can live several Gyrs, and even remain 
strong over a Hubble time. This can solve the problem of 
maintenance of this perturbation in some isolated galaxies. 



6 CONCLUSIONS 

By treating the lopsided distribution in an exponential 
galactic disc as a global feature, and using the softened self- 
gravity of the perturbation, we show that the resulting lop- 
sided mode is long lived to several Gyr (essentially increasing 
the winding up time scale), and it also exhibits a radially in- 
creasing amplitude as observed. The emerging pattern speed 
of the lopsided mode is extremely slow being a factor of ~ 
10 smaller than the local free precession rate and with no 
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particularly preferred sense of precession. The behaviour of 
the pattern speed is confirmed in the numerical simulation. 

Our semi-analytical global treatment shows that the 
lopsidedness appears to be a purely oscillatory normal mode 
outside the continuum in smaller sized disc, while it emerges 
as an instability in comparatively larger sized disc with a 
very slowly precessing pattern. The e-folding growth time 
scale for a typical disc like that of Milky way is ~ 1.2 Gyr. 

Our numerical analysis show that the disc is never un- 
stable to the lopsided modes without the indirect term aris- 
ing due to the lopsided perturbation itself. Even though 
there is a subtle interplay between the indirect term due 
to the lopsided perturbation and the collective effect of the 
softened self-gravity, it appears that the indirect term plays 
a crucial role in preparing the disc susceptible to lopsided 
instability. 

Numerical simulations of an isolated purely stellar com- 
ponent, an exponential disc without any stellar bulge or 
dark matter halo, has confirmed the theoretical predictions, 
and also show long-lived m = 1 global modes. The criteria 
for stability are different than for local WKB modes, and 
gravity softening stabilises essentially the local modes, but 
the global ones are not much affected by the softening. The 
m = 1 modes superpose to the bar mode, that appear in 
all runs studied here. The buckling instability forming the 
peanut shape combines to the lopsidedness to produce disc 
tilts. 

Further work will include spheroidal collisionless com- 
ponents, and also the gas component in the disc, in order 
to compare more realistically these results with the obser- 
vations. 
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